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Fate decision processes of T lymphocytes are crucial for health and disease. Whether aT 
lymphocyte is activated, divides, gets anergic, or initiates apoptosis depends on extracellu- 
lar triggers and intracellular signaling. Free cytosolic calcium dynamics plays an important 
role in this context. The relative contributions of store-derived calcium entry and calcium 
entry from extracellular space toT lymphocyte activation are still a matter of debate. Here 
we develop a quantitative mathematical model of T lymphocyte calcium dynamics in order 
to establish a tool which allows to disentangle cause-effect relationships between ion fluxes 
and observed calcium time courses. The model is based on single transmembrane protein 
characteristics which have been determined in independent experiments. This reduces the 
number of unknown parameters in the model to a minimum and ensures the predictive 
power of the model. Simulation results are subsequently used for an analysis of whole cell 
calcium dynamics measured under various experimental conditions. The model accounts 
for a variety of these conditions, which supports the suitability of the modeling approach. 
The simulation results suggest a model in which calcium dynamics dominantly relies on 
the opening of channels in calcium stores while calcium entry through calcium-release 
activated channels (CRAC) is more associated with the maintenance of theT lymphocyte 
calcium levels and prevents the cell from calcium depletion. Our findings indicate that 
CRAC guarantees a long-term stable calcium level which is required for cell survival and 
sustained calcium enhancement. 
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1. INTRODUCTION 

The dynamics of the free cytosolic calcium concentration upon 
stimulation of T lymphocytes (TCs) is crucial for TC activation 
and fate decision processes. While it is clear that calcium rises 
upon stimulation of the TC receptor (TCR), the calcium pattern 
associated with different fates of TCs has not been deciphered ( 1- 
4). However, it is likely that the calcium signal is correlated with 
the later fate of the activated TC, i.e., anergy, division, acquisition 
of the regulatory phenotype or apoptosis (3, 5, 6). Dysregulation 
in TC calcium signaling has been linked to inflammatory and 
autoimmune diseases as well as to allograft rejection (2). A rele- 
vant player in calcium dynamics is the calcium-release activated 
channel (CRAC) which is located in the plasma-membrane and 
activated by calcium depletion in the intracellular calcium stores 
like the endoplasmic reticulum (ER). As ion-gating transmem- 
brane proteins in the plasma-membrane (PM) are possible targets 
of drug applications in the context of various clinical settings (2, 
7, 8), insight into the specific calcium dynamics is essential for an 
efficient control of TC behavior. 

The relative contribution of ER-derived calcium entry versus 
CRAC to the calcium signal following TC stimulation is a mat- 
ter of ongoing debate (9-12). While a considerable number of 
scientists argue for CRAC being the major player of TC calcium 



dynamics (12, 13), others argue for a dominant role of calcium- 
induced calcium -release (CICR) (10, 14, 15) or for a dominant 
role of second messenger-induced calcium-release from ER (16, 
17). All three contributions are required for a functional TC cal- 
cium signal, however, the sequence of the contributions might 
be essential. Second messenger-induced activity appears as a very 
early signal (18), which might act as triggering event for CICR 
and subsequent CRAC activation (2, 19-21). Quantitative analysis 
of the components of calcium signaling during TC activation is 
essential for the development of strategies for an efficient control 
of TC responses. A mathematical analysis of the calcium dynamics 
in a model including calcium stores and CRAC may shed light on 
the relation and relevance of both calcium sources (12, 22). This 
is the major motivation for the present work. 

T lymphocytes are non-excitable cells in the sense that they 
do not exhibit bursts like pancreatic beta-cells, spikes like neu- 
rons, or comparable fast dynamics (23-25) even though single 
cell measurements detected calcium oscillations (10, 19, 26). The 
non-excitability of TCs might have prevented a larger interest of 
mathematical modelers in lymphocyte calcium dynamics. The few 
existing models (22, 27-29) mostly focused on modeling of CRAC- 
channel dynamics or a special part of the store-operated-calcium- 
entry signaling pathway (27) and its contribution to intracellular 
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FIGURE 1 | Scheme of the transmembrane proteins included in the 
mathematical model. The elements of aTC considered in the 
mathematical model are shown with particular emphasis on 
calcium-conducting transmembrane proteins. The outer border of the cell 
represents the PM. In the PM CRAC and PMCA are located which induce 
calcium influx and extrusion, respectively. The intracellular organelle around 
the nucleus (yellow) depicts the ER.The membrane of the ER contains 
IP3R and SERCA (sarco/endoplasmic reticulum calcium-ATPase) which 
control the exchange of calcium between the ER and the cytosol. 



calcium dynamics. Also the dependence of ORAI1 assembly to 
a tetrameric CRAC on calcium oscillations was considered (28). 
In one approach a spatial resolution of CRAC currents and of 
the calcium dynamics in TCs was considered in the context of 
immunological synapse formation (22). Two different models for 
inositol 1,4,5-tris-phosphate-receptor (IP3R) activity were com- 
pared and it was shown that they differ in their impact on TC 
calcium dynamics (27), a result that will be used in the present 
model as well. The plasma-membrane calcium-ATPase (PMCA) 
was modeled in Jurkat TCs and a reversible modulation of PMCA 
activity was postulated (12), which is a further topic addressed in 
this investigation. 

The aforementioned theoretical studies on TC calcium dynam- 
ics were all based on whole cell current models. In the context 
of excitable cells we have shown that it is possible to derive the 
whole cell currents from the single transmembrane properties 
(30). To achieve this goal, specific quantitative measurements of 
protein activation, inactivation, dependencies and conductance 
were used. The big advantage of this approach is that most para- 
meters of the model are determined by independent experiments, 
which increases the predictive power of the model. The present 
paper applies this strategy to calcium dynamics of TCs. The model 
also includes the dynamics of transmembrane channel expression 
kinetics in order to represent CRAC recruitment upon calcium 
store depletion. The mathematical model was validated using 
dynamic calcium data measured under specific experimental con- 
ditions. The validated model allowed to reassess the relative role 
of store-derived and CRAC-mediated calcium entry on a quanti- 
tative basis. A new role of the CRAC-channel is postulated, which 
is associated with maintenance of TC calcium levels rather than 
TC activation. 

2. MATERIALS AND METHODS 

The modeling framework is presented in this section. Three com- 
partments, extracellular space (ES), cytosol, and ER are distin- 
guished, each being represented by ordinary differential equations. 
The nucleus is only included as an object which reduces intra- 
cellular space (Section 2.4). The compartments are encased by 
the PM and the membrane of the ER. Both membranes contain 
transmembrane proteins (Figure 1), which allow for a flow of 
ions from one compartment to the other. The surface densities 
and the properties of these proteins in terms of conductance and 
control parameters determine the resting and activated states of 
the TC. While the protein properties are derived from measured 
single protein data, which are independent of the whole cell cal- 
cium experiments used for validation, the protein densities are 
in parts derived from steady state conditions (Section 2.9) and 
in parts determined by fitting to whole cell calcium dynamics 
(Section 2.10). 

In the following, the equations for calcium dynamics, calcium- 
buffering, and for the second messenger D-myo-inositol- 1,4,5- 
tris-phosphate (IP3) are introduced. The details of the com- 
partment sizes and the surface between the compartments are 
explained in Section 2.4. Particular emphasis is put on the geom- 
etry of the lurkat TC, the specific cell line, which is used in the 
subsequently described experiments. The single transmembrane 
protein characteristics will be described and the corresponding 



mathematical models introduced. Wherever possible we imple- 
mented data from experiments performed with Jurkat TCs. Finally, 
all model pieces are merged to the proposed model of TC calcium 
dynamics in Section 2.10. This includes the determination of the 
remaining unknown parameters. 

2.1. CALCIUM DYNAMICS AND BUFFERING 

Two major players determine exchange of calcium through the PM 
(Figure 1): PMCAs actively transport calcium from the cytosol 
into ES (12, 31), while CRAC-channels allow for a passive elec- 
trochemical current of calcium ions into the cell (21, 32, 33). The 
density of active CRAC-channels Pcrac in the PM is increased 
in dependence on ER calcium (Cer) depletion (34). The free 
cytosolic calcium concentration (C) is further affected by modu- 
lations of the calcium flow between cytosol and ER (Figure 1). 
Sarco/endoplasmic reticulum calcium-ATPase (SERCA) trans- 
ports calcium ions from the cytosol into the ER (10, 35) and by 
this maintains a chemical gradient of calcium from the ER to the 
cytosol. Conversely, calcium can passively leave the ER into the 
cytosol when IP3R channels open in dependence on calcium and 
the second messenger IP3 (1, 36, 37) associated with calcium- 
induced-calcium-release (CICR) (38). Furthermore other second 
messengers like cyclic ADP ribose (cADPR) (16) and nicotinic 
acid adenine dinucleotide phosphate (NAADP) (39) were found to 
influence calcium dynamics. Their effect is mediated by activation 
of the ryanodine receptor (RyR) which leads to calcium-release 
from intracellular stores (17, 40-43). In order to avoid an over- 
parametrization of whole cell calcium curves, which the present 
model focusses on, we restrict ourselves to the dynamics of IP3. 
The inclusion of cADPR and NAADP and their effect on RyR 
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requires a more detailed data basis and should be addressed with 
a more complex model in the future. 

2. 1. 1. Calcium in the cytosol 

The four sources and sinks of free cytosolic calcium C are described 
as 



dC 



-1 



(£PPMCAJ~PMCA + fPCRAC^CRAC 



dt 2t a F(l + B c ) 

+ ?ERCPSERCAJ~SERCA + ?ERCPlP3R^IP3R) > (1) 

where zc a = 2 is the valence of calcium ions, F the Faraday 
constant, £ and £erc geometrical surface to volume ratios for PM 
[equation (16)] and ER membrane [equation (17)], respectively. 
Px is the surface density and 7x the single transmembrane protein 
current, which are defined in the subsequent sections. By conven- 
tion, positive ions that enter the cytosol are represented by negative 
electrical currents. 5c represents the cytosolic calcium-buffer in 
the rapid buffer approximation 



Bc = 



b 0 K h 



(c + K h y 



(2) 



where bo is the total buffer concentration, and K\, the calcium- 
buffer dissociation constant. The fraction of free calcium in the 
cytosol then reads 



fc = 



1 + 



C + K h 



(3) 



The main buffer within the cytosol is calmodulin (CaM) with 
4 calcium binding sites per CaM. There is a diversity of mea- 
sured CaM concentrations depending on cell type and organ 
(44-46). A realistic average value is 25 /xM of CaM, corresponding 
to bo = 100 /xM of calcium binding sites. The dissociation con- 
stant K\, was determined by the required fraction of free calcium 
of fc ~ 0.1% in non-excitable cells (47), leading to K\, = 0.1 /U-M. 

2. 1.2. Calcium in the ER 

The dynamics of the calcium concentration in the ER Cer is 
described by an equation analogous to equation (1) 



dCER ?ER (PSERCAJSERCA + PIP3RJIP3R) 



dt 



ZCaF (1 + BCER) 



(4) 



with a different geometrical factor §er [equation (18)], and a 
different calcium-buffer Bqer defined by 



5c,ER = 



^ER,0^ER,b 



(Cer + KER,b) 

The fraction of free calcium in the ER reads 



[ Cer + -Ker,i,] 



-l 



(5) 



(6) 



The resting ER calcium level is Cer,o = 400 /xM (2) which holds 
true for Jurkat TCs considered here (34). 



In the ER calcium is buffered by calsequestrin, with three high 
and three low-affinity binding sites (48), and by calreticulin, with 
two distinct domains, one with high-affinity (K = 0.0l mM) but 
low capacity (0.6-1 mol Ca 2+ /mol protein), and one with low- 
affinity (K = 2 mM) but high capacity (18 mol Ca 2+ /mol protein) 
(49). As calreticulin binds more than 50% of the luminal ER cal- 
cium (50) only this buffer is considered here. In pancreatic acinar 
cells it was estimated that 20-times more calcium would be free 
in the ER compared to the cytosol (47), suggesting /c,er ~ 0.02. 
This is achieved by using foER,o = 30mM and K^n^ = 0.1 mM, 
which corresponds to an intermediate dissociation constant of 
both calcium binding domains. 

2.2. IP3 DYNAMICS 

IP3 (P) is generated in a TCR- and calcium-dependent way 
described by 



dP 
~dt 



= p P H(C,C v ,n P )T(t)-yeP , 



(7) 



where flp is the production and yp the degradation rate. The 
Hill-function is defined as 



H(X,K, n) 



X" + K n 



(8) 



where K is the concentration of X at which the Hill-function 
reaches its half value, and n the Hill-coefficient which determines 
the steepness of the Hill-function. 

The degradation rate in equation (7) is determined by steady 
state conditions for IP3 in equation (35). The production rate is 
the tonic production rate and is modulated by increased calcium 
with the Hill-function in equation (7), leading to a positive feed- 
back loop between calcium and IP3. jSp is fitted as described in 
Section 2.10 and mainly influences the speed of the early calcium 
response after TC stimulation. 

The production is further modified by the time-dependent 
input function T{t) representing the degree of TCR stimulation 
of the cell. T = 1 is assumed in the resting state. 

The resting concentration of IP3 Pq is identified as critical 
parameter. It strongly determines the responsivity of the cell via 
activation of IP3R (see Section 2.6). It was fitted as described in 
Section 2.10. 

2.3. MEMBRANE AND REVERSAL POTENTIALS 

The resting membrane potential is set to V = — 60 mV (51-53). 
It is assumed that the membrane potential is not changed by the 
calcium currents (V = Vq) and that the electrical current corre- 
sponding to calcium fluxes in or out of the cell is equilibrated by 
other ions. 

Further it is assumed that Ver,o = Vo = —60 mV, thus, the ER 
and the cytosol are electrically equilibrated (54). ER calcium efflux 
may lead to small fluctuations (55) which are neglected. Thus, 
Ver = V is assumed at all times. 

In this approximation, the reversal potentials depend on the 
chemical gradient only. The Nernst-equation is used to calcu- 
late the reversal potential during dynamical changes of calcium 
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concentrations: 



_ RT , 
V c = — - In 



(Qxt\ 
~c)~ 



Vc,ER 



rt ^ 



AV C 



In 1-^1 - A Vqer 



(9) 



where _R = 8.315 J/(K mol) the Rydberg (molar) gas constant, 
T = 310 K, and F the Faraday constant. 

For many channels, the real I-V-relationship is not linear as 
assumed for the currents 7x in equations (23) and (28). There- 
fore, the reversal potential is corrected for the CRAC-channel by a 
shift A Vc = 78 mV in order to achieve the correct linear extrap- 
olation of the I-V-relation of CRAC-channels with a zero around 
Vc ^ 50 mv [(34) Figure 1, (33) Figure 2], This approximation is 
only valid for depolarization below V = 50 mV. 

The reversal potential for ER calcium Vqer is treated in com- 
plete analogy to the cytosolic case which leads to a correction of 
A Vqer- As the value is not known it is derived using the fitting 
routine in Section 2.10. 

2.4. TC GEOMETRY 

Most measurements on TC calcium dynamics are performed in 
Jurkat TCs which are small but still larger than normal human 
blood derived TCs. In an approach based on ordinary differential 
equations the effect of an ion current onto the concentration of 
the ion in the cytosol or ER is not spatially resolved. While local 
calcium entry can induce transient high concentrations of cal- 
cium (22) the comparably small cytosolic volume of TCs justifies 
this approach for the description of whole cell calcium dynam- 
ics because local inhomogeneities will quickly equilibrate. In the 
model this is reflected as change of the average concentration. 
How an ion current changes the average calcium concentration 
depends on the geometry of the cell. In the dynamic equations for 
the ion concentrations the current 7x through an individual ion- 
conducting protein X is multiplied by the surface density px- Thus 
the concentration change is derived from a current surface density. 
The latter has to be translated to the actual change in concentration 
by a surface to volume ratio, which is considered here. 

Given a cell radius _R C ell> the cell volume V ce ii and cell surface 
A ce u are known as well. However, the volume relevant to changes in 
concentration is not the cell volume V ce u but the cytosolic volume 
V cyt which can be approximated as 



V< 



cyt 



Vcell - Ver - V nucleus 



(10) 



using the volumes of ER and nucleus. This is important because 
the nucleus, with a radius of 



^nucleus — /R^cell 



(ID 



is substantially reducing the resulting cytosolic volume, /r ~ 0.8 is 
assumed for human TCs (56), and /r & 0.25 for Jurkat TCs. The 
volume of the ER is expressed as a fraction of the total cell volume 



with/v ~ 0.1 (57). However, electron micrographs of TCs suggest 
that /v~0.01 (25) which is used here. Taking this together, the 
cytosolic volume becomes 



V cyt = V cell (l-/ V -/ R 3 ) 



(13) 



The surface of the ER is also needed in order to translate 
the current surface densities calculated on the ER surface into 
concentration changes in cytosol and ER. While the TC itself is 
approximated as a sphere, the ER is absolutely non-spherical. The 
exact surface of the ER is difficult to be measured and accordingly 
approximated as 



A E R = JaAer = 4:71 f A 



3 Ver 
4tt 



2/3 



(14) 



where Aer is the surface of a spherical ER with volume Ver [deter- 
mined in equation (12)], and /a is the fold increase of the ER 
surface with respect to the surface of a spherical ER. /a = 30 was 
roughly estimated from the folding degree of the ER in elec- 
tron micrographs of TCs (25). Note that only the product of 
/A with fy^ 3 enters the model, such that both parameters are 
redundant and were only separated because of their physiological 
meaning. 

The size of human blood TCs can be estimated starting from 
the capacity of C m = 0.028 pF//xm 2 (58, 59) and using the whole 
cell capacitance of C ce u = 2 pF [(60), p. 606]. C ce ii = 1.7 pF was 
found in Fomina et al. (14). Using C C eii/A ce ii = C m this yields a 
radius 



Rcett = 



Qell 
4jrC m 



(15) 



and the resulting _R ce ii ~ 2.4 /xm corresponds to A ce ii = 72.4 /im 2 . 
The experiments described below were performed with Jurkat 
TCs and the same authors determined the cell volume to 
V ce ll = 2pl (12). This determines the values J? ce u = 8/im and 
A ce n = 804.2 /im 2 used in the present simulations. 

Given the cell radius £ C ell> the fractions jv and/R, as well as the 
factor /a, the surface to volume ratios required in equations (1) 
and (4) can be calculated by 



§ERC : 



Acell 

V cyt 

A E r 
V cyt 



. A ER 



(16) 



(17) 



(18) 



with 



Ver = jv V ce ii 



(12) 



Acell = 4nR 



cell 



(19) 



Frontiers in Immunology |T Cell Biology 



September 2013 | Volume 4 | Article 277 | 4 



Schmeitz et al. 



Modeling calcium dynamics of T lymphocytes 



Vc Y t=-7tR 3 cell (l-fv-fi) (20) 
V E R=^^ c 3 ell (21) 

/ ~ \ 2/3 

A ER = 4jtf A ( ^ j . (22) 
2.5. CRAC-CHANNEL 

The open CRAC-channel current is determined by the electro- 
chemical gradient 

Jcrac = gCRAC (V - Vq) ■ (23) 

This approach closely follows the model of Martin et al. (22). 
The validity of the Ohm's law approximation is only guaranteed 
within limited ranges of membrane potentials. 

2.5.1. Single channel conductance 

The single channel CRAC conductance was found to be extremely 
small in the order of gcRAC = 2 fs (61). 

2.5.2. CRAC recruitment 

The density of active CRAC-channels, estimated by the steady state 
CRAC-channel current, is a dynamic function of the ER-calcium 
concentration [(34) Figure 1C] described by 

^PCRAC _ PCRAC - PCRAC . . 

77 — ' (24j 

at TCRAC 

where 

PCRAC = PcRAC + (PCRAC ~ PcRAc) 

(1 — H (Cer, Ccrac> »crac)) , (25) 

with Ccrac = 169 /xM and «crac = 4.2. To our knowledge, this is 
the first time that the surface density of active CRAC-channels in 
the PM is modeled as a dynamic quantity. 

When estimating the same quantity from the degree of 
STIM1 -redistribution, a rather similar relation is found with 
Ccrac = 187 /xM and mcrac = 3.8 [(34), Figure 2]. The unifor- 
mity of both curves supports the view that CRAC-channels are 
recruited and open in response to ER calcium depletion (34). 
It can be deduced from the similarity of both curves that the 
opening dynamics is substantially faster than the redistribution 
of STIM1. As no opening dynamics of the CRAC-channel is 
included in the model, the dynamics of the current itself and not 
of STIM1 -redistribution is used. 

2.5.3. CRAC density 

In equation (25) Pcrac are tne u PP er an d lower limits of possi- 
ble active CRAC densities. The resting value Pcrac,o is not known 
from experiment and is determined by parameter fitting to calcium 
dynamics upon TCR stimulation (Section 2.10). 

The density of CRAC-channels upon activation with PHA 
increased about 9-fold (33) which constraints Pqraq- ^ 10-fold 
increase has been reported for the whole cell CRAC current in 



response to stepwise reduced Cer [Figure 1C in Luik et al. (34)]. 
These findings translate into the condition 

Pcrac = /cracPcraqo • (26) 

where Pqrac was determined by parameter fitting within the 
experimental boundaries in Section 2.10. A value for Pq^q 
is not known and is determined by the steady state condition 
equation (36). 

2.5.4. CRAC time scales 

The time scale of CRAC recruitment can be estimated from the 
rising time of calcium curves which provides an upper bound of 
tcrac < 100 s for the activation time. It is likely that this time is 
associated with CRAC recruitment rather than with CRAC open- 
ing because opening time scales are typically much shorter. The 
time scale of CRAC recruitment is set to tcrac = 5 s. Larger values 
could also be used as the fit was insensitive to tcrac- Inactivation 
of CRAC-channels is difficult to be assessed (62). As the time 
scale of inactivation is in the order of 1000 s (14) and thus longer 
than the typical experimental durations used here, the present 
model ignores CRAC inactivation and assumes that the reduction 
of active CRAC-channels is a secondary effect of Cer recovery. 

2.6. IP3R IN THE ER 

The ER releases its calcium content if activated by IP3 and cy tosolic 
calcium (63, 64). The release of calcium from intracellular stores 
is based on the opening dynamics of RyR and IP3R in the mem- 
brane of the ER. TCs express both, RyR and IP3R and even different 
subtypes of them. 

IP3Rs have binding sites for IP3 and calcium and exhibit com- 
plex forms of cooperativity (65). For the present purpose the 
heuristic description of IP3R activation and inhibition is sufficient. 
The characteristic feature of the IP3R conductance is a calcium- 
dependent log-bell-shaped opening probability curve (63) which 
has been measured for ER vesicles from canine cerebellum and 
further reviewed in Foskett et al. (38). The open probability curve 
was best fitted by the product of an activating and an inhibiting 
Hill-function, both with Hill- coefficient 2 (63). 

2.6.1. Open probability 

TCR signaling leads to the generation of IP3, the ligand of the IP3R, 
which modulates the open probability of IP3R and is described as 
the product of an activation term gn>3R and an inactivation term 
^IP3R- The properties of single channel openings were quantita- 
tively determined in Xenopus laevis oocytes (37, 38) and we assume 
that the single channel properties are transferable to TCs. The mea- 
sured dynamics are well described by the previously published 
Mak-McBride-Foskett model (37, 38). 

gIP3R = gIP3R,maxH (C, ClP3R jact , «IP3R,act) 
felP3R = H (Cip 3 R ; inh, C, «IP3R,inh) 
ClP3R,inh = Qp3R,inh^ (P, ^IP3R,C> «IP3R,c) ■ (27) 

The according parameters were obtained from a data fit 

(37) to be gIP3R,max = 0.81, ClP3R,act = 0.21/i,M, ttTp 3R;act = 1 .9, 
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FIGURE 2 | Activation and inactivation of IP3R in dependence on calcium 
and IP3. Reproduction of the experimental data (37) with the 
Mak-McBride-Foskett model equation (27). (A) IP3R opening probability 



giP 3R h| P3R in dependence on free cytosolic calcium C for different 
IP3-concentrations P. (B)The same IP3R opening probability in dependence 
on IP3-concentration P for different free cytosolic calcium concentrations C. 



«iP3R,inh = 3.9, Cip 3Riin h = 52 fj.M, «iP3R,c = 4, and the IPs- 
concentration of half-activation Pip3r,c = 0.050 fiM. 

The dependencies of the model equation (27) on calcium and 
IP3 are depicted in Figure 2 and correctly reproduce the measure- 
ments in Mak et al. (37) suggesting that the modulating effect of 
IP3 is mediated by IP3R inactivation ( 38) . At low calcium this effect 
is hardly visible and IP3R activation remains unaffected by changes 
in IP3 for resting concentrations beyond 50 nM. The dynamic IP3 
range is between 100 nM and 1 /xM (see (66), Table 1), a regime 
of IP3 at which the IP3R-typel exhibits saturation (27, 67). We 
do not aim at resolving whether the resting IP3 is lower in TCs 
or whether the IP3R characteristics are different in TCs, such that 
the DeYoung-Keizer model should be employed instead (27, 65). 
It is assumed that the resting concentration of IP3 is in the range 
of 5-10 nM which ensures that an increased IP3-concentration 
has an impact on the IP3R opening probability as depicted in 
Figure 2. 

2.6.2. IP3R calcium current 

The steady state activation function (Figure 2) can be used to 
define the calcium current through the IP3R which follows the 
electrochemical gradient between the cytosol and the ER 



I\P3R = gIP3RgIP3R^IP3R(V — Ver — Vqer)) (28) 

with jip3r = 0.064pS. Note that the conductance differs between 
tissues (38). V — Ver is the potential difference between ER and 
cytosol, and Vqer is the ER-reversal potential for calcium, cal- 
culated from the Nernst-equation equation (9). We assume an 
electrical equilibrated relation of cytosol and ER such that V= Ver 
holds true. 

2.6.3. (Injactivation dynamics 

The activation and inactivation factors g and h are treated dynam- 
ically and approach equation (27) in steady state, while the 



adaptation of C]P3R ; i n h in equation (27) is treated in quasi steady 
state: 

dglP3R _ glP3R,maxH(C, ClP3R jact , «IP3R,act) — gIP3R 

dt TIP3R 
<#»IP3R _ #(Qp3R,inh> C, «IP3R,inh) ~ ^IP3R 

dt # IP3R 

2.6.4. Activation time 

Activation time scales can be determined from Mak and Foskett 

(68) , Figure 5, and are in the range of less than 5 and 20 ms for 
depolarizations to 20 and 40 mV. Two exponentials were needed to 
fit the opening frequency. Using rat hepatocytes the activation and 
inactivation time scales were found to depend on the IP3-levels 

(69) : activation varied between 100 and 500 ms for 10 /xM and 
400 nM of IP3, respectively (see Figure 1 therein). The time delay 
reported in Marchant and Taylor (69) is consistent with the IP3- 
dependent time delay of channel opening of 1 s > tip3r > 100 ms 
using basophilic leukemia cells from rats (70). As the model 
focusses on calcium dynamics on the scale of minutes, a constant 
TIP3R = 100 ms is assumed. 

2.6.5. Inactivation time 

Onset of inactivation happens in less than 2 min (68). A slow and 
a fast current were distinguished (69). The fast current inactivates 
on a time scale of 200-450 ms, the slow one between 1 and 6 s (see 
Figure 2 in the same publication). The authors attribute the fast 
time scale to inactivation of IP3R and the slow one to the deple- 
tion of the calcium content in the ER. Accordingly, only the fast 
time scales are relevant for the single IP3R, and #ip3r = 300 ms is 
assumed. 

2.6.6. Calmodulin dependence 

It was found that the calcium-release from ER is reduced for 
high concentrations of the calcium-buffer calmodulin (71). Such 
a dependence is neglected in the present model. 
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2.6 J. IP3R density 

The IP3R density on Jurkat TCs is not known and is determined 
using steady state condition equation (34). 

2.7. PLASMA-MEMBRANE CALCIUM-ATPase 

Plasma-membrane calcium-ATPase is an ATP-driven calcium 
pump which extrudes calcium from the cell to the ES. It was char- 
acterized in TCs (12). In a first attempt the dependence on the 
ATP concentration is ignored and assumed to be large enough in 
order to make the pump work optimally. In this case the activity 
is mainly dependent on the calcium concentration in the cell. A 
suitable modeling approach is 



with 



JpMCA = JpMCAgPMCA 



%PMCA H (C, CpMCA> «PMCa) — gPMCA 



dt 



^PMCA 



(30) 



(31) 



The current Jpmca is positive as it carries calcium out of the 
cell. The Hill- coefficient was determined to be mpmca = 2 (72). 

2.7.1. Turn-overrate 

The turn-over rate of PMCA is approximately 30 Hz (73) which 
corresponds to an activity rate of /c a = 0.03/ms which is also 
used in Sherman et al. (74). This turn-over rate can be trans- 
lated into an electrical current by using that every pumping 
event corresponds to the flow of 2 electrical charges which yields 
5mca = ZCa^ a = 60- 1.6- 10~ 19 C/s « 10~ 17 A = 10~ 5 pA. 

2.7.2. Calcium-dependent activation 

Typical values for the half- activation calcium concentration are 
Cpmca = 0.1 /xM (at 540 nM calmodulin, see e.g., (75), Figure 3). 
The isoforms 2a and 2b exhibit CpMCA2ab < 0.1 /xM. The predom- 
inant isoform of PMCA in Jurkat TCs is 4b [(76), Figure 6B]. 
CpMCA4b ^ 0.2 /xM was found in Jurkat TCs [Figure 2 in Caride et 
al. (76)] and is used here. 

2.7.3. Calmodulin-dependent activation 

Note that the values of half- activation also depend on the 
calmodulin concentration (75-77). For calmodulin concentra- 
tions above 1 /xM (which is exclusively the case in all present 
simulations) full activation of all isoforms is ensured (75, 78). 
Hence, the dependence on calmodulin is weak in this regime and 
is neglected. 

2.7.4. Delay of activation 

Binding of calcium to PMCA is a comparably fast process with a 
rate constant >3 per second (79). However, the activity of PMCA 
is delayed in some isoforms including the isoform 4b (76) which 
is relevant for TCs. The rate constant of PMCA activation upon 
stimulation with 500 nM of free calcium was in the range of 0.02 
per second (76), which leads to a delay of PMCA activation in the 
range of tpmca = 50 s in equation (31). This delay was associated 
with a calmodulin and calcium-dependent activation (76). How- 
ever, as the exact mechanism is not known this delay is modeled 
in equation (31) on a phenomenological basis. 



2.7.5. PMCA density 

For TCs no precise value of the protein density is known and the 
value is determined by the steady state condition equation (33). 

2.8. SARC0/END0PLASMIC RETICULUM CALCIUM-ATPase 

The calcium level in ER is kept high with the help of SERCA cal- 
cium pumps. The activity of SERCA is assumed to rapidly adapt 
to the present calcium concentration in the cytosol and can be 
described by a Hill-function. 



Jserca = jsercaH(C, Cserca> "serca)- 



(32) 



In every turn-over cycle two calcium ions are transported per 
ATP (80). There are different subtypes of SERCA whose properties 
differ. In Jurkat TCs as well as in human tonsil lymphocytes the 
dominant isoform is SERCA2b (81). 

2.8.1. Turn-overrate 

The turn-over rates of most isoforms are in the range 
of /c=10Hz (i.e., j"sERCA = asERCAZCa^=6- 10~ 6 pA). For 
SERCA2b k 2 b = 5 Hz was reported (82) which implies the value 
Jserca = 3 • 10~ 6 pA used in the model. 

2.8.2. Calcium-dependent activation 

For the SERCA isoforms 1, 2a, and 2b the half-activation 
Cserca = 0.4 /xM and the Hill-coefficient mserca = 2 are a good 
approximation [(82) Figure 4]. The half-activation of SERCA3 
is around 1/xM with the same Hill-coefficient (82). SERCA2b is 
an isoform active at relatively low calcium concentrations (82). 
Specifically, in Jurkat TCs as well as in human tonsil lympho- 
cytes, the dominant isoform SERCA2b was characterized with 
«SERCA2b = 2.0 and C S ERCA2b = 0.25 /xM (81), which is used here. 

2.8.3. SERCA density 

Even though the expression of SERCA was shown to be mod- 
ulated upon activation (83), the expression density of SERCA 
within ER is not known and is determined by parameter fitting 
in Section 2.10. 

2.9. STEADY STATE DETERMINES PROTEIN DENSITIES 

The resting state of the TC is determined by setting the dynamics 
in equations 1, 4, 7, and 24 to zero. Accordingly, the equations for 
C, C E r, P, and p C RAC, read: 



Ppmca = — 



PIP3R 



yp : 



PC RAC,oJcRAC,0 
JPMCA.O 
_ PSERCA^SERCA.O 
-flP3R,0 

fipH(C 0 , C P , n P ) 



Po 



(33) 
(34) 
(35) 



^CRAC — 



PCRAQO 



H (Cer,o, Ccrao «crac) 
[l - /crac (l - H (Cer,o, Ccrac. «crac))] . (36) 

where 7x,o denotes the currents with all quantities X in the resting 
configuration Xo. The parameters of the model are summarized 
in Table 1. 
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Table 1 | List of parameters, values, and references. 



Parameter 


Meaning 


Value 


AQI (%) 


Comments and reference 






f?cell 


TC radius 


8/xm 


4.0 


Fixed from JurkatTC (12) 


fn 


Nucleus radius fraction of R ce ii 


0.25 


0.1 


Fixed, Section 2.4 


fy 


ER volume fraction 


0.01 


2.8 


Fixed (25) 


fA 


ER surface, fold of spherical 


30 


3.5 


Fixed, Section 2.4 


Cm 


Membrane capacitance 


28fF//im 2 


0.0 


Fixed (58, 59) 




T 


Temperature 


310K 




Fixed, used in equation (9) 


V 0 


Resting membrane potential 


-60 mV 




Fixed (51, 52, 53) 


Ver,o 


Resting ER potential 


-60 mV 




Fixed = V o (54) 


Co 


Resting calcium 


0.1 fiM 




Fixed (12) 


Cer.o 


Resting ER calcium 


0.4mM 


4.0 


Fixed in range 0.1-0.8 mM (2) 


Cex\ 


Extracellular calcium 


2mM 


0.0 


Fixed by experimental settings 


AV C 


Reversal potential shift 


78 mV 


0.3 


Fixed (33, 34) 


A^CER 


ER-reversal potential shift 


63 mV 


176 


Variable 




bo 


Cytosolic calcium-buffer 


100 /xM 


3.6 


Fixed (44, 45, 46) 


K b 


Buffer dissociation constant 


0.1 ,tM 


0.2 


Fixed by f c = 0.1% in equation (3) 


£>ER,0 


ER calcium-buffer 


30 mM 


4.3 


Fixed by f c ,ER = 20f c (47) 


^ER,b 


ER buffer dissociation constant 


0.1 mM 


1.2 


Fixed (49) 


Po 


Resting IP3 


8.7 nM 


318 


Variable 


P> 


IP3 production rate 


0.6nM/s 


6.5 


Variable 


yp 


IP3 degradation rate 


0.01149/s 




Fixed by steady state equation (35) 


c P 


Calcium of half IP3 production 


0.5/iM 


19.4 


Fixed (85) 


n P 


IP3 production Hill-coefficient 


1 


21.2 


Fixed 


TRANSMEMBRANE PROTEIN DENSITIES 


PIP3R 


ER-IP3R density 


11.35//xm 2 




Fixed by steady state equation (34) 


PSERCA 


ER-SERCA density 


700///m 2 


3.7 


Variable 


PPMCA 


PMCA density 


68.57//xm 2 




Fixed by steady state equation (33) 


PCRACO 


Resting active CRAC density 


0.6/ynm 2 


1.0 


Variable 


^CRAC 


Max active CRAC density 


3.9//im 2 


6.5 


Variable, range from Luik et al. (34) 


^CRAC 


Min active CRAC density 


0.5115//nm 2 




Fixed by steady state equation (36) 



Fixed parameters were determined either by direct measurement, by indirect constraints, or using steady state conditions. Variable parameters were subject to the 
fitting algorithm described in Section 2. 10. AQI measures the sensitivity of Ql in equation (37) for changed parameter values: each parameter is increased by 10% 
and the percentage of the change in Ql is provided. Values below 0.05% are given as 0.0%. 



2.10. NUMERICAL SOLUTION AND PARAMETER FITTING 

The model defined by the equations (1, 2, 4, 5, 7, 8, 16-25, 
27-32) was implemented as C ++ -code and solved using a self- 
written 4th-order Runge-Kutta algorithm with adaptive stepsize 
control. 

As not all parameters could be determined by steady state 
conditions or by experimental constraints, Figure 1A in Bautista 
et al. (12) was used to determine the remaining free para- 
meters. We used a two-step fitting procedure: at first, the 
differential evolution algorithm defined in Storn and Price 
(84) was incorporated into the C ++ -code of the model on 
the basis of all parameters in Table 1 that were not deter- 
mined by steady state conditions. The parameters were varied 
within hard-coded boundaries dictated by experimental con- 
straints (when available). The quality of the fit to the calcium 



data in Bautista et al. (12) was measured as the mean square 
deviation 



QI = 

N, 



\ ,=i 



(X, - Eif 



(37) 



with X[ and E{ representing the simulation and experimental 
values, respectively. In a second step, the first approximative 
fit was subject to a sensitivity analysis in which each parame- 
ter was varied by 10% while monitoring the effect on Ql. The 
three unknown protein densities pserca, Pcraqo, and Pcrac 
two sensitive IP3 related parameters Po and /Jp, as well as 
the very sensitive parameter AVqer were used for fine-tuning 
the initial parameter fit with the same differential evolution 
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algorithm. These fit parameters are marked variable in Table 1. 
The final fit reached QI = 0.100197 with N = 22. The sensi- 
tivity analysis was repeated for the final fit and the impact 
of each parameter on QI in equation (37) is provided in 
Table 1. 

All subsequently described simulations are started with the cell 
in steady state as defined by the hard-coded equations (33-36). 
Starting from these initial conditions, the respective stimulation 
protocols are applied as described in the results section. 

3. RESULTS 

In the methods section single protein characteristics were sum- 
marized and specific mathematical models capturing their main 
properties were proposed or cited. The models for the single pro- 
teins were combined to a whole cell model and the unknown 
parameters were determined using steady state conditions or by 
data fitting as described in Section 2. 10. In this section, we replicate 
specific experimental setups described in the literature in silico and 
analyze the calcium dynamics from the perspective of the model. 

3.1. TCR STIMULATION 

The introduced TC-model is used to investigate the experiments 
of Bautista et al. (12) in Jurkat TCs. TC activation by stimula- 
tion of TCR induces an intracellular rise of second messengers like 
cADPR, NAADP, and IP3. In the present model this rise is col- 
lectively reflected in equation (7) for IP3. The IP3 signal activates 
IP3R and by this induces a calcium-release from the ER. The pos- 
itive feedback loop of CICR leads to even more calcium-release 
from the ER, which in turn reduces the ER calcium concentration 
Cer. CRAC is activated in a CER-dependent manner, as repre- 
sented by equation (24). The rising cytosolic calcium is cleared by 
PMCA and the ER is refilled by SERCA, both being ATP-dependent 
processes. 

Using 2 mM of external calcium a cell in resting state was acti- 
vated with OKT3 via TCR [see Figure 1A in Bautista et al. (12)]. 
This induced a calcium peak to more than 1 ^iM within 50-100 s, 
which subsequently relaxed to a plateau level of ~0.7 /xM over the 
following 200 s. This behavior is well reproduced by the model 
(Figure 3A) and was used to determine the unknown parameters 
(Table 1). The peak height relies on both, on a proper activation 
of IP3R and CRAC. The choice of Po turned out to be rather 
important in order to guarantee a proper activation of IP3R. The 
maximum activated CRAC density was essential for the height 
of the plateau. A value of /crac = 6-5 was found to correctly 
reproduce the plateau height as measured in Bautista et al. (12). 

Bautista et al. reported that the delay of PMCA activation is 
responsible for the calcium overshoot ( 12). Accordingly, we inves- 
tigated whether this conclusion is supported by the model. For 
that purpose we reduced tj>mca from 50 s (76) to 1 millisecond 
in the otherwise unaltered simulation. This modification does not 
touch the steady state configuration, such that all other parameters 
remained unchanged. The model still generated an overshoot of 
calcium, however, with a reduced amplitude (Figure 3A, blue dot- 
ted line). The height of the plateau as well as the relaxation time 
from peak to plateau remained unchanged. Thus, the simulation 
supports an influence of the PMCA delay onto the overshoot, but 
it turned out to not be essential for its existence. 



This surprising result led to the question whether a differ- 
ent delay in the model could explain the calcium overshoot. All 
delays were tested and the only delay with impact on the over- 
shoot was IP3R inactivation, i.e., the parameter 6>ip3r (Figure 3A, 
green dashed line) . In conclusion, the model suggests that intrinsic 
properties of the IP3R and not the delay of PMCA activation are 
responsible for the calcium overshoot upon TCR stimulation. 

Next, the model is used for an analysis of the currents which 
are associated with the different phases of the calcium dynamics 
(Figure 3B). The IP3-current (red full line) is clearly the largest 
and also persists beyond the overshoot of calcium due to SERCA 
activity (red dotted line). The calcium peak induces a PMCA cur- 
rent (black dotted line) that drives the calcium out of the cell (black 
dashed line). Thus, the CRAC current (black full line), induced by 
the loss of calcium in the ER, does not even induce a net flow 
of calcium into the cell (black dashed line) but just prevents the 
cell from running out of calcium. This model result suggests that 
the role of CRAC is the stabilization of the cell rather than its 
activation, which is mostly mediated by calcium from the ER. 

The single protein currents (Figure 3C) show that in the model 
IP3R currents react much more dynamic than CRAC currents. 
The increased cytosolic calcium even reduces the CRAC currents 
on the single channel level. However, the reduced ER calcium 
level induces a strong increase in the active CRAC-channel density 
(Figure 3D). Thus, the increased whole cell CRAC current seen 
in Figure 3B is not a result of single channel responses but of a 
changed channel density. 

3.2. TCR STIMULATION WITH ZERO EXTERNAL CALCIUM 

Zero external calcium experiments aim at suppressing CRAC cur- 
rents in order to investigate ER calcium currents in response to 
TCR stimulation. TCR stimulation of Jurkat TCs hold at zero 
external calcium results in an intracellular calcium peak after 50- 
100 s with reduced amplitude in comparison to stimulation with 
normal external calcium [Figure 2A in Bautista et al. (12)]. The 
increased calcium is cleared below baseline level within 100-200 s. 

In the model, the Nernst-equation prohibits the use of zero 
external calcium conditions. At very low calcium concentrations 
the Nernst-equation loses its validity and the reversal potential 
diverges. Therefore, external calcium is set to the concentration 
C* xt at which the CRAC current vanishes in resting state: 

Qxt = Co exp j — j . (38) 

This mimicks the suppression of CRAC currents as intended 
in the experiment. The in silico result is shown in Figure 4. The 
CRAC current vanishes at resting state (Figure 4B, black full line). 
This is approximately true during the whole experiment, such that 
the method for mimicking the zero external calcium experiment 
appears appropriate. 

The calcium peak (Figure 4A) is lower than the one found in 
Figure 3 which is qualitatively consistent with the experimental 
result (12). Furthermore, the time scales of calcium rise and clear- 
ance are perfectly matched between simulation and experiment. 
Even the overshoot of clearance below the baseline calcium level 
is fully reproduced. However, quantitatively, the peak is higher in 
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FIGURE 3 | Calcium dynamics in response toTCR stimulation Simulation 
of the experiment in Figure 1 A in Bautista et al. (12). (A) The cell is in steady 
state until f= 10s when stimulation is started by setting 7~(f)= 1.6 in equation 
(7). Stimulation is kept constant throughout the simulation. Reference is the 
simulation used as a basis for all other simulations in the paper and is 
compared to the experimental result [black triangles read off Figure 1A in 



Bautista et al. (12)]. For fast PMCA a value of 1 ms was used for t PM ca. For 
slow IP3R inactivation a value of 300 s was used for fl| P3H . (B) Whole cell 
currents show an initial release of calcium from the ER which is followed by a 
CRAC inward current. Negative currents are calcium currents into the cytosol. 
(C) Single transmembrane protein currents. (D) Dynamic response of active 
CRAC-channel density. 




FIGURE 4 | Calcium dynamics in response toTCR stimulation at zero 
external calcium. Simulation of the experiment in Figure 2A in Bautista et al. 
(12). (A) The cell is in steady state until f= 10s, when stimulation is started by 
setting 7(f) = 1.6 in equation (7) and external calcium is set to equation (38). 



Stimulation and external calcium are kept constant throughout the simulation. 
(B) Whole cell currents show the almost complete inhibition of CRAC 
currents. The calcium peak (left) is generated by ER calcium only. Negative 
currents are calcium currents into the cytosol. 
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theory than in experiment. As the in silico peak is exclusively gen- 
erated by the ER, one might hypothesize that the ER is too big or 
its calcium content is too high. As these parameters were already 
chosen comparably low (Table 1), it is more likely that the lack of 
external calcium influences the stimulation of the cell. The peak 
size can be reduced to the measured amplitude by reducing the 
stimulation T from 1.6 to 1.25 (not shown). 

3.3. BLOCK OF SERCA AT ZERO EXTERNAL CALCIUM 

Thapsigargin (TG) is frequently used to block SERCA activity as 
it prevents calcium uptake by the ER and leads to a continuous 
reduction of ER calcium. As low ER calcium recruits and activates 
CRAC-channels, this would lead to a strong influx of extracel- 
lular calcium. In order to prevent this according ER depletion 
and CRAC activation experiments are performed in zero calcium 
medium. This strategy was used in Jurkat TCs to generate a cell 
state in which the ER is mostly void of calcium and CRAC-channels 
are recruited and activated to a maximum (12, 13, 25, 86). It was 
reported that this procedure leads to a transient calcium peak 
of about 0.5 /iM after more than 100 s which is slowly cleared 
and reaches calcium levels below the resting level [Figure 1A in 
Quintana et al. (86)]. 

Having established a strategy [equation (38)] for mimicking a 
medium with zero calcium, a SERCA block is performed in silico 
by setting /serca = 0 at f = 10 s. No TCR stimulation was applied. 
The measured dynamics are well reproduced without any fur- 
ther parameter fitting (Figure 5A, black full line). Furthermore, 
the intended depletion of ER is achieved (Figure 5A, red dashed 
line): a continuous reduction of ER calcium is observed. Note that 
also IP3 exhibits some dynamics (Figure 5A, blue dotted line) 
which further accelerates ER calcium loss by activation of IP3R. 
As expected, the reduced ER calcium leads to the recruitment of 
active CRAC-channels (Figure 5B). 

3.4. BLOCK OF PMCA IN A TG TREATED TC 

As the TG-mediated SERCA block works in silico (Figure 5), the 
role of PMCA in the clearance of cytosolic calcium is investigated. 



Following Figure 6 in Bautista et al. (12) TG was applied in a zero 
calcium medium as in Figure 5. Then a pulse of 2 mM external 
calcium was applied for 50 s. This induces a steep rise in calcium 
which is also steeply cleared again. Together with Lanthan (La 3+ ), 
a PMCA, and CRAC inhibitor (2), the clearance of such a calcium 
peak was substantially slower ( 1 2 ) . A block of PMCA alone by car- 
boxyeosin led to an only weakly modified clearance time, without 
a return to baseline levels within 300 s. 

The same protocol is applied in the model (Figure 6). As before, 
ER calcium is depleted (Figure 6A, red dashed line) and the active 
CRAC density is increased correspondingly (as in Figure 5B). The 
initial free cytosolic calcium peak (Figure 6A, black full line) is 
the same as in Figure 5A. Upon increasing external calcium to 
2 mM for 50 s at t = 300 s, a strong CRAC current is induced 
(Figure 6B, black full line), which steeply increases cytosolic cal- 
cium (Figure 6A, black full line). This calcium is also quickly 
cleared upon return to the mimicked zero calcium medium. Cal- 
cium clearance is dominated by the PMCA current (Figure 6B, 
black dotted line). However, in the model a small CRAC current is 
observed supporting extrusion of calcium out of the cell in the case 
of zero external calcium concentration. Upon permanent block of 
PMCA and repetition of the transient stimulation by external cal- 
cium, it is this backward CRAC current that clears calcium from 
the cytosol. The time course of the clearance is slower than without 
PMCA block and the calcium baseline is not reached after 350 s 
(Figure 6B, black full line). This is a similar behavior as in the cor- 
responding experiment with carboxyeosin [Figure 6D in Bautista 
et al. ( 1 2 ) ] . However, it is not known whether the real CRAC allows 
for such inverse current under zero calcium conditions. The return 
of calcium to the baseline might also be supported by uptake of cal- 
cium by other organelles like mitochondria, which is not covered 
by the present model. 

In the case of PMCA- and CRAC-block with La 3+ a 
rather slow calcium clearance is observed in experiments 
which apply the same stimulation protocol [Figure 6C in 
Bautista et al. (12)]. In silico no clearance is observed at 




FIGURE 5 | Calcium dynamics in response toTG SERCA block at 
zero external calcium. Simulation of the experiment in Figure 1 A in 
Bautista et al. (86). (A) The cell is in steady state until f= 10s, when 
SERCA block is applied and external calcium is set to equation (38). 



Block and external calcium are kept constant throughout the 
simulation. P (blue dotted line), C (black full line), and C ER (red dashed 
line) are shown. Note the factor 1000 applied to P and C. (B)The time 
course of active CRAC density pcrac- 
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FIGURE 6 | Calcium dynamics in response to SERCA and PMCA block 
at zero external calcium. Simulation of the experiment in Figures 6C,D of 
Bautista et al. (12). (A) The cell is in steady state until f= 10s, when 
SERCA block is applied and external calcium is set to equation (38). At 
f= 300s 2 mM external calcium are applied for 50 s followed by a switch 
back to equation (38). The procedure is repeated at r=600s. In addition, 
/pmca is set to zero for the remaining simulation time. SERCA block is kept 
throughout the simulation. P (blue dotted line), C (black full line), and C ER 



(red dashed line) are shown. Note the factor 1000 applied to P and C. 
(B)The CRAC density was increased during the first 300s by ER calcium 
depletion with TG at zero external calcium [(A) red dashed line]. The whole 
cell currents show a sudden CRAC current (black full line) upon restoration 
of external calcium to 2 mM, which is cleared by PMCA activity (black 
dotted line) after return to zero external calcium. When PMCA is blocked in 
addition (at f = 600 s), the calcium clearance is slower (black full line). 
Negative currents are calcium currents into the cytosol. 
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FIGURE 7 I Block of PMCA in untreated TCs. TheTCR stimulation protocol 
in Figure 3 is repeated. At the time of stimulation (f = 10s), PMCA activity 
is blocked preventing calcium efflux from the cell. Cytosolic calcium (full 
black line, left axis) and the response of the CRAC density (dotted red line, 
right axis) are shown. The axis for the CRAC density was scaled as in 
Figure 3D for better comparison. 



all. SERCA is blocked and positive 7iP3R-currents are not 
allowed in the model, such that an uptake of calcium into 
the ER is excluded. A full block of PMCA and CRAC also 
excludes any expulsion of calcium out of the cell. Thus, the 
model suggests that the slow clearance of cytosolic calcium 
is either due to incomplete block of PMCA or to leakage 
currents. 

3.5. BLOCK OF PMCA IN AN UNTREATED TC 

The model suggests that the role of CRAC for TC activation 
is mainly the maintenance of the integrity of the TCs during 
stimulation in the sense that it prevents the activated TC from 
running out of calcium. If we block PMCA in silico at the time 
of TCR stimulation (protocol as in Figure 3) in an otherwise 
untreated TC, the TC would be prevented of losing calcium. 
According to our interpretation of the role of CRAC we would 
expect that CRAC activity is strongly reduced in comparison to 
Figure 3. 

T lymphocytes receptor stimulation of PMCA-blocked TCs 
is predicted to induce a strong increase of cytosolic calcium 
(Figure 7, full black line). Thereby, the steady state CRAC cur- 
rent is not increased (as in Figure 3B, full black line) but reduced 
(not shown). However, the CRAC current is not reduced to zero 
such that the total block of PMCA infers a persistent net influx of 
calcium into the cell and, thus, to a persistent increase of cytoso- 
lic calcium. This would ultimately destroy a real TC. As SERCA 
activity is normal in this in silico experiment, calcium in the ER 
is only transiently reduced (not shown), leading to a transient 
and weak increase in the CRAC density (Figure 7, dotted red 
line). The amplitude is less than twofold instead of sixfold in 
Figure 3D. 



This result supports the interpretation of the role of CRAC- 
channels in silico and may be tested in experiment. It further shows, 
that the calcium level in the ER may be used as an indicator for 
the overall calcium status of TCs. 

4. DISCUSSION 

Within the presented investigation whole cell calcium dynam- 
ics were derived from models of single transmembrane 
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ion-conducting proteins. The model successfully described 
a number of experimental settings and captured impor- 
tant characteristics of calcium dynamics upon TCR stim- 
ulation. In particular, the role of store-derived calcium- 
release versus CRAC activation was well represented. We con- 
clude, that we have generated a modeling framework suit- 
able for the quantitative analysis of calcium dynamics of 
TCs. 

The value of using measured single transmembrane protein 
characteristics is twofold: at first, it substantially reduces the num- 
ber of free parameters in an otherwise very complex model. All 
measured single protein properties were just implemented and 
not altered in the fine-tuning of the model. The reduced number 
of free parameters increases the predictive power of the model. 
Secondly, it can be considered as a multi-scale approach which 
links single proteins to whole cell behavior. This allows the analy- 
sis of the dynamics on the single protein level and its implications 
onto the cellular properties. The drawback of this approach is 
that not all parameters could be derived from TCs such that we 
had to assume that single transmembrane properties are univer- 
sal, which is a correct assumption in many cases (30). Within this 
framework, cell-specific properties are controlled by the protein 
densities and the cell-specific properties, both determining the 
activity range of the respective proteins. However, the predictive 
power of the model would benefit from corresponding TC related 
data. 

The model of Jurkat TC calcium dynamics as supported by the 
mathematical model starts from a TCR-derived increase in second 
messengers like IP3. This ultimately activates IP3R currents and 
induces an initial calcium current from the ER into the cytosol 
which triggers CICR. However, calcium uptake via SERCA activ- 
ity equilibrates the calcium loss in the ER, which leads to a zero 
net flux on long-term - despite ongoing TCR stimulation. The 
ER calcium loss, according to the model, was the major contribu- 
tion to the free cytosolic calcium peak. However, as PMCA activity 
leads to an overall loss of calcium in the cell, a compensation 
mechanism is required for a sustained elevation of free cytosolic 
calcium. The model suggests that the CRAC activity essentially 
contributes to this compensation mechanism. The initiation of 
CRAC currents by the depletion of ER calcium levels is in line 
with this interpretation. The model predicts, that a standard stim- 
ulation of TCs together with a block of PMCA activity would 
lead to a strong calcium rise together with a minor and tran- 
sient increase of the CRAC density and on long-term to a reduced 
CRAC current (see Figure 7). This prediction may be tested in 
experiment in order to validate this interpretation of the role of 
CRAC-channels. 

The emerging hypothesis that the role of CRAC is the stabi- 
lization of the TC calcium level, needs to be further strengthened 
by more detailed modeling work. In particular, early events after 
stimulation like NAADP generation (87) and subsequent RyR cal- 
cium currents from the ER are essential for the early calcium 
rise (11, 18) and will have to be included in the mathematical 
model for a proper time-resolved coverage of calcium dynam- 
ics. The need for additional mechanisms is also underpinned by 



the strong sensitivity of the model behavior to changes in the 
IP3 resting concentration Pq (Table 1). In the present simula- 
tion the long-term calcium plateau height mostly relies on the 
maximum possible CRAC activation by store-operated calcium 
depletion. For the long-lasting calcium rise, cADPR was proven 
relevant (16, 40) and has to be considered in the context of the 
present hypothesis on the role of CRAC. Also the transferability 
to human blood derived TCs, which exhibit a different geometry, 
has to be assessed. 

As the PMCA block experiment at zero external calcium 
has shown, it might be important to include leakage currents 
into the model. However, it should be noted that the pos- 
tulated role of delayed PMCA activity for the calcium over- 
shoot after TCR stimulation (12) could only partially be con- 
firmed by the mathematical model. With fast activation of 
PMCA, an overshoot was still observed in our model, and 
the overshoot could only be suppressed by lack of IP3R 
inactivation. It should be noted that this result might rely 
on the Mak-McBride-Foskett model (37), which was used 
for IP3R dynamics and which exhibits inactivation at rather 
low IP3. The result might differ if the TC calcium dynam- 
ics would be based on the DeYoung-Keizer model for IP3R 
activity (65). 

The proposed model has limitations in its range of applic- 
ability. For example, the usage of the Nernst-equation for the 
chemical gradient and Ohm's law for the current-voltage rela- 
tionship of ion-conducting pores is justified only in narrow limits. 
Experiments with zero external calcium drive the model to the 
very limits of this range of applicability which was circumvented 
here using a phenomenological approximation. The model maybe 
reformulated in terms of the Fokker-Planck-equation in order to 
describe ion transport through the pores in more detail. Further- 
more, it is known that calcium entrance points lead to spatially 
inhomogeneous calcium dynamics (22) which are not covered 
by the present space-averaged model. The value of the present 
approach lies in the surprising result that quantitative characteris- 
tics of single transmembrane proteins are sufficient to determine 
the cell behavior in the framework of an ordinary-differential- 
equation based model. The model has proven its predictive power, 
as it was fitted to data of one experiment in Figure 3 and could 
be used to predict and explain further data of calcium dynam- 
ics generated under other experimental conditions in Figures 4 
and 6. It is planned to elaborate the potential and the limits 
of the model predictions by application to further experimental 
settings. 
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